function frfrefch = sinefrf02
%   sinefrf function computes a set of frequency response functions based on swept sine 
%   time history data.
%
%Author:    Donald J. Hershfeld
%   ManTech International Corporation
%Created:   October 8, 2002
%Released:
%Purpose: The sinefrf function requests the user to select a channel of data to be used
%         as the frequency reference and a channel to be used as the frequency response 
%         function reference.The frequency reference channel is analyzed using a Hilbert
%         transform to determine the rate of phase angle change with respect to time. 
%         The user is then prompted to select break points on the phase angle plot.
%         Logarithmically swept cosine and sine references are then fitted between the
%         break points. 
%
%         By multiplying each time history by the cosine and sine references, exactly half
%         of the energy at the reference frequency is transposed to zero frequency. The other half of
%         the energy is transposed to twice the reference frequency. A low pass filter is
%         used to create a tracking filter. Two times the output of the filter is equal to
%         the peak amplitude of the time history of the reference frequency. The cosine part
%         is the in-phase component and the sine part is the quadrature component. Each 
%         frequency response function is stored as a complex array in the Analysis structure.
%
%         The frequency response function reference channel is analyzed first. All of the
%         channels are then analyzed and divided by the reference channel to obtain each FRF.
%
%Function Parameters:
%   At this time, the sinefrf function does not require input parameters. However,
%   the sinefrf function could be changed at a later time to require a freqrefchan and a frfrefchan.
%
%   freqrefchan:  
%       Channel number of data defined by Project.Test_Item.Run.Channel(ChanNum)to be used to
%       estimate the instantaneous rate of change of phase angle (reference frequency).
%
%   frfrefchan:   
%       Channel number of data defined by Project.Test_Item.Run.Channel(ChanNum)to be used as
%       the reference for all freqeuncy response functions.
%
global Project;
refch = sineswpdialog;
frfrefch = sinerefdialog;
TimeRef=Project.Test_Item.Run.Time_Ref;
% Get the freqeuncy reference data.
% Modified 14-Apr-2004 by CM: Convert single to double precision
x=double(Project.Test_Item.Run.Channel(refch).Time_His);
h=hilbert(x);  %  Compute the Hilbert transform.
PhaseAngle=atan2(imag(h),real(h));  %  Calculate the instantaneous phase angle.

%  Calculate the rate of change of phase angle.
SampleRate=Project.Test_Item.Run.Sample_Rate;
freq = SampleRate * diff(unwrap(PhaseAngle))/(2*pi);
[B,A]=butter(4,8/Project.Test_Item.Run.Sample_Rate);
filtfreq=filter(B,A,freq);
% Added 03-May-2007 CMcLeod to eliminate values < 0 in semilogy()
filtfreqP = filtfreq;
for ii = 1:size(filtfreq,2)
    if (filtfreq(ii) < 0)
        filtfreqP(ii) = 0;
    else
        filtfreqP(ii) = filtfreq(ii);
    end
end
semilogy(TimeRef(2:length(TimeRef)),filtfreqP),grid on;  %   Display the instantaneous frequency.
% set(gca,'ylim',[1 250]);
set(gca,'ylim',[1 3000]);
if (Project.Test_Item.Run.Time_Ref(length(TimeRef)) < 180)
    set(gca,'XTick',0:5:Project.Test_Item.Run.Time_Ref(length(TimeRef)));
else
    set(gca,'XTick',0:10:Project.Test_Item.Run.Time_Ref(length(TimeRef)));
end
xlabel('Time - seconds');ylabel('Frequency - Hz.');
%fprintf('Select breakpoints on figure.'); 

% uiwait(msgbox(mycellarray,' ','modal'));
% uiwait(msgbox(mycellarray,'Breakpoints','modal'));
%mycellarray={'Select breakpoints on figure.'};
%Hf_BP = msgbox(mycellarray,'Breakpoints','modal');
% get(Hf_BP)
%uiwait(Hf_BP);

% [start stop]=start_stop_dialog;
status = 0;
while status == 0  %  Endpts not acceptable
    [xi,yi]=ginput(2);
	x = [min(xi) max(xi)];
	haxes = get(gca,'Children');
	xd = get(haxes,'xdata');
	yd = get(haxes,'ydata');
	xlow = find(x(1)>xd);
	xhigh = find(x(2)<xd);
	starttime = xd(size(xlow,2));
	startfreq = yd(size(xlow,2));
	stoptime = xd(xhigh(1));
	stopfreq = yd(xhigh(1));
	[status] = ConfirmEndpts(starttime,startfreq,stoptime,stopfreq);
end  % Endpts acceptable
first=round(SampleRate*starttime);last=round(SampleRate*stoptime);% used to be start and stop

RefFreq = filtfreq(first:100:last);     %  Decimate the reference frequency.
%  Calculate the referrence cosine and sine signals.
cx = 2 * cos(PhaseAngle);
sx = 2 * sin(PhaseAngle);
% disp( 10/Project.Test_Item.Run.Sample_Rate );      % 200 percent?
[B,A]=butter(2,8/Project.Test_Item.Run.Sample_Rate);   %[B,A]=butter(4,.005);  %  The low-pass filter.

% Compute sine reference channel first
ich=frfrefch;
refr=filtfilt(B,A, double(Project.Test_Item.Run.Channel(ich).Time_His(first:last)) .* cx(first:last));
refi=filtfilt(B,A, double(Project.Test_Item.Run.Channel(ich).Time_His(first:last)) .* sx(first:last));
hwait = waitbar(0,'Processing Sine Data: Please wait...');
Project.Test_Item.Run.Num_Chan;

% Determine Which Analysis index to save Data to.
try
 	n = length(Project.Test_Item.Run.Channel(1).Analysis);
 	n = n + 1;
catch
	n = 1;
end
   
for ich=1:Project.Test_Item.Run.Num_Chan
%     str1=Project.Test_Item.Run.Channel(ich).Trans_Loc;
%     str2=Project.Test_Item.Run.Channel(ich).Coordinate;
%     str3=Project.Test_Item.Run.Channel(ich).Direction;
%     str4=Project.Test_Item.Run.Channel(ich).Orientation;
    yr=filtfilt(B,A, double(Project.Test_Item.Run.Channel(ich).Time_His (first:last)) .* cx(first:last));
    yi=filtfilt(B,A, double(Project.Test_Item.Run.Channel(ich).Time_His (first:last)) .* sx(first:last));
    h=(yr+i*yi)./(refr+i*refi);
    h_mod = h(1:100:length(h));
    Project.Test_Item.Run.Channel(ich).Analysis(n).Type=5; 
    Project.Test_Item.Run.Channel(ich).Analysis(n).Abs=RefFreq;
    Project.Test_Item.Run.Channel(ich).Analysis(n).Ord=h_mod;
	Project.Test_Item.Run.Channel(ich).Analysis(n).FPoint = first;
	Project.Test_Item.Run.Channel(ich).Analysis(n).LPoint = last;
	Project.Test_Item.Run.Channel(ich).Analysis(n).RefChan = frfrefch;
    waitbar(ich/Project.Test_Item.Run.Num_Chan);
end

close(hwait);